home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / shgeqz.z / shgeqz
Text File  |  1996-03-14  |  10KB  |  265 lines
  1.  
  2.  
  3.  
  4. SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      SHGEQZ - implement a single-/double-shift version of the QZ method for
  10.      finding the generalized eigenvalues  w(j)=(ALPHAR(j) +
  11.      i*ALPHAI(j))/BETAR(j) of the equation   det( A - w(i) B ) = 0  In
  12.      addition, the pair A,B may be reduced to generalized Schur form
  13.  
  14. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  15.      SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB,
  16.                         ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
  17.                         INFO )
  18.  
  19.          CHARACTER      COMPQ, COMPZ, JOB
  20.  
  21.          INTEGER        IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
  22.  
  23.          REAL           A( LDA, * ), ALPHAI( * ), ALPHAR( * ), B( LDB, * ),
  24.                         BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
  25.  
  26. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  27.      SHGEQZ implements a single-/double-shift version of the QZ method for
  28.      finding the generalized eigenvalues B is upper triangular, and A is block
  29.      upper triangular, where the diagonal blocks are either 1-by-1 or 2-by-2,
  30.      the 2-by-2 blocks having complex generalized eigenvalues (see the
  31.      description of the argument JOB.)
  32.  
  33.      If JOB='S', then the pair (A,B) is simultaneously reduced to Schur form
  34.      by applying one orthogonal tranformation (usually called Q) on the left
  35.      and another (usually called Z) on the right.  The 2-by-2 upper-triangular
  36.      diagonal blocks of B corresponding to 2-by-2 blocks of A will be reduced
  37.      to positive diagonal matrices.  (I.e., if A(j+1,j) is non-zero, then
  38.      B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be positive.)
  39.  
  40.      If JOB='E', then at each iteration, the same transformations are
  41.      computed, but they are only applied to those parts of A and B which are
  42.      needed to compute ALPHAR, ALPHAI, and BETAR.
  43.  
  44.      If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal
  45.      transformations used to reduce (A,B) are accumulated into the arrays Q
  46.      and Z s.t.:
  47.  
  48.           Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
  49.           Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*
  50.  
  51.      Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
  52.           Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
  53.           pp. 241--256.
  54.  
  55.  
  56. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  57.      JOB     (input) CHARACTER*1
  58.              = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will not
  59.              necessarily be put into generalized Schur form.  = 'S': put A and
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  71.  
  72.  
  73.  
  74.              B into generalized Schur form, as well as computing ALPHAR,
  75.              ALPHAI, and BETA.
  76.  
  77.      COMPQ   (input) CHARACTER*1
  78.              = 'N': do not modify Q.
  79.              = 'V': multiply the array Q on the right by the transpose of the
  80.              orthogonal tranformation that is applied to the left side of A
  81.              and B to reduce them to Schur form.  = 'I': like COMPQ='V',
  82.              except that Q will be initialized to the identity first.
  83.  
  84.      COMPZ   (input) CHARACTER*1
  85.              = 'N': do not modify Z.
  86.              = 'V': multiply the array Z on the right by the orthogonal
  87.              tranformation that is applied to the right side of A and B to
  88.              reduce them to Schur form.  = 'I': like COMPZ='V', except that Z
  89.              will be initialized to the identity first.
  90.  
  91.      N       (input) INTEGER
  92.              The order of the matrices A, B, Q, and Z.  N >= 0.
  93.  
  94.      ILO     (input) INTEGER
  95.              IHI     (input) INTEGER It is assumed that A is already upper
  96.              triangular in rows and columns 1:ILO-1 and IHI+1:N.  1 <= ILO <=
  97.              IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
  98.  
  99.      A       (input/output) REAL array, dimension (LDA, N)
  100.              On entry, the N-by-N upper Hessenberg matrix A.  Elements below
  101.              the subdiagonal must be zero.  If JOB='S', then on exit A and B
  102.              will have been simultaneously reduced to generalized Schur form.
  103.              If JOB='E', then on exit A will have been destroyed.  The
  104.              diagonal blocks will be correct, but the off-diagonal portion
  105.              will be meaningless.
  106.  
  107.      LDA     (input) INTEGER
  108.              The leading dimension of the array A.  LDA >= max( 1, N ).
  109.  
  110.      B       (input/output) REAL array, dimension (LDB, N)
  111.              On entry, the N-by-N upper triangular matrix B.  Elements below
  112.              the diagonal must be zero.  2-by-2 blocks in B corresponding to
  113.              2-by-2 blocks in A will be reduced to positive diagonal form.
  114.              (I.e., if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and
  115.              B(j,j) and B(j+1,j+1) will be positive.)  If JOB='S', then on
  116.              exit A and B will have been simultaneously reduced to Schur form.
  117.              If JOB='E', then on exit B will have been destroyed.  Elements
  118.              corresponding to diagonal blocks of A will be correct, but the
  119.              off-diagonal portion will be meaningless.
  120.  
  121.      LDB     (input) INTEGER
  122.              The leading dimension of the array B.  LDB >= max( 1, N ).
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      ALPHAR  (output) REAL array, dimension (N)
  141.              ALPHAR(1:N) will be set to real parts of the diagonal elements of
  142.              A that would result from reducing A and B to Schur form and then
  143.              further reducing them both to triangular form using unitary
  144.              transformations s.t. the diagonal of B was non-negative real.
  145.              Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0),
  146.              then ALPHAR(j)=A(j,j).  Note that the (real or complex) values
  147.              (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized
  148.              eigenvalues of the matrix pencil A - wB.
  149.  
  150.      ALPHAI  (output) REAL array, dimension (N)
  151.              ALPHAI(1:N) will be set to imaginary parts of the diagonal
  152.              elements of A that would result from reducing A and B to Schur
  153.              form and then further reducing them both to triangular form using
  154.              unitary transformations s.t. the diagonal of B was non-negative
  155.              real.  Thus, if A(j,j) is in a 1-by-1 block (i.e.,
  156.              A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.  Note that the (real or
  157.              complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are
  158.              the generalized eigenvalues of the matrix pencil A - wB.
  159.  
  160.      BETA    (output) REAL array, dimension (N)
  161.              BETA(1:N) will be set to the (real) diagonal elements of B that
  162.              would result from reducing A and B to Schur form and then further
  163.              reducing them both to triangular form using unitary
  164.              transformations s.t. the diagonal of B was non-negative real.
  165.              Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0),
  166.              then BETA(j)=B(j,j).  Note that the (real or complex) values
  167.              (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized
  168.              eigenvalues of the matrix pencil A - wB.  (Note that BETA(1:N)
  169.              will always be non-negative, and no BETAI is necessary.)
  170.  
  171.      Q       (input/output) REAL array, dimension (LDQ, N)
  172.              If COMPQ='N', then Q will not be referenced.  If COMPQ='V' or
  173.              'I', then the transpose of the orthogonal transformations which
  174.              are applied to A and B on the left will be applied to the array Q
  175.              on the right.
  176.  
  177.      LDQ     (input) INTEGER
  178.              The leading dimension of the array Q.  LDQ >= 1.  If COMPQ='V' or
  179.              'I', then LDQ >= N.
  180.  
  181.      Z       (input/output) REAL array, dimension (LDZ, N)
  182.              If COMPZ='N', then Z will not be referenced.  If COMPZ='V' or
  183.              'I', then the orthogonal transformations which are applied to A
  184.              and B on the right will be applied to the array Z on the right.
  185.  
  186.      LDZ     (input) INTEGER
  187.              The leading dimension of the array Z.  LDZ >= 1.  If COMPZ='V' or
  188.              'I', then LDZ >= N.
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          SSSSHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  203.  
  204.  
  205.  
  206.      WORK    (workspace/output) REAL array, dimension (LWORK)
  207.              On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
  208.  
  209.      LWORK   (input) INTEGER
  210.              The dimension of the array WORK.  LWORK >= max(1,N).
  211.  
  212.      INFO    (output) INTEGER
  213.              = 0: successful exit
  214.              < 0: if INFO = -i, the i-th argument had an illegal value
  215.              = 1,...,N: the QZ iteration did not converge.  (A,B) is not in
  216.              Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i=INFO+1,...,N
  217.              should be correct.  = N+1,...,2*N: the shift calculation failed.
  218.              (A,B) is not in Schur form, but ALPHAR(i), ALPHAI(i), and
  219.              BETA(i), i=INFO-N+1,...,N should be correct.  > 2*N:     various
  220.              "impossible" errors.
  221.  
  222. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  223.      Iteration counters:
  224.  
  225.      JITER  -- counts iterations.
  226.      IITER  -- counts iterations run since ILAST was last
  227.                changed.  This is therefore reset only when a 1-by-1 or
  228.                2-by-2 block deflates off the bottom.
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.